function synapse()
%
% calculate and make sample data for a binomial release process.
%
% define parameters:
% m - quantal content
%

expt = 'original';
% or
%expt = 'treated';
% original values
k=[0:15]; % range for binomial histogram
factor = 5;
N = 10000;
switch(expt)
    case 'original'
        m=8; % quantal content
        mEPSC = 10;
        tau = 1; % rate for alpha function

    case 'treated'
        % part d values:
        m = 4;
        mEPSC = 7;
        tau = 5;
    case 'test'
        m = 1;
        mEPSC = 1;
        tau = 1;
    otherwise
        return;
end;


noise  = 2.2; % noise, in "pA" for current trace

nk = binomial(m, k, N);
newfigure('synapse');

im = k * factor;

subplot('Position', [0.1, 0.35, 0.8, 0.5]);
bar(im, nk);
xlabel(gca, 'EPSC (pA)');
ylabel(gca, '# of observations');

mc = sum(nk .* im)/N;  % get mean value (should be close to mEPSC)
for i = 1:length(nk)
    fprintf(1, '%6d\t%6d\n', im(i), floor(nk(i)));
end;

pk = nk./N; %expressed as probability

sepsc2 = sqrt(sum(pk.*im.^2)-mc^2); % an equatin I found that is wrong
sepsc = sqrt(sum(pk.*(im - mc).^2)); % the correct equation 
%(Weisstein, Eric W. "Variance." From MathWorld--A Wolfram Web Resource.
%http://mathworld.wolfram.com/Variance.html
cv = sepsc/mc;
mest = 1/(cv^2);
cv2 = sepsc2/mc;
mest2 = 1/(cv2^2);

% from method of failures:
mfail = log(N/nk(1));

% print out the values
fprintf(1, 'Mean of epsc: %7.3f  Mean of histogram: %7.3f\n', mEPSC*factor, mc);
fprintf(1, 'True m = %7.2f   N = %d\n', m, N);
fprintf(1, 'CV m   = %7.2f   (sd = %7.3f)\n', mest, sepsc);
fprintf(1, 'CV2 m  = %7.2f   (sd = %7.3f)\n', mest2, sepsc2);
fprintf(1, 'MFailm = %7.2f   \n', mfail);
%
% Now plot some EPSCs for comparison...
%
sr = 0.1;
fs = 1000/sr;
lpf = 1000;
lpfn = 2500;

a = randexpo(0.05, 25);
a = cumsum(a); % turn into cumulative times
tx = [0:sr:max(a)+100];

ai = floor(a/sr); % find when events occur
zt=zeros(length(tx), 1);
zt(ai) = mEPSC;
zt = zt .* (randn(length(zt), 1) * 0.15 + 1);
af = alpha(tau, [0:sr:100]);
epsc = zeros(1, length(zt));
epsc = conv(zt, af); % convolution to get alpha function
[bn,an] = butter(8,lpfn/(fs/2));
Hn = dfilt.df2t(bn,an);          % Direct-form II transposed structure
epscn = filter(Hn,randn(length(epsc), 1) * noise); % filter the noise too
epsc = epsc + epscn; % add broadband noise
% then filter the data
[b,a] = butter(8,lpf/(fs/2));
Hd = dfilt.df2t(b,a);          % Direct-form II transposed structure
epsc = filter(Hd,epsc);



% the result is actually longer than zt...
tx = [0:sr:(length(epsc)-1)*sr];


subplot('Position', [0.1, 0.1, 0.8, 0.15]);
plot(tx, epsc);
ylabel(gca, 'Im (pA)');
xlabel(gca, 'time (msec)');
set(gca, 'XLim', [0,600]);




function [nk] = binomial(m, k, N)
nk = N*m.^k .* exp(-m)./factorial(k);





function [a] = alpha(tau, x)
a =  x ./ tau .* exp(1 - x ./ tau);
